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Abstract. We introduce a method for determining the functional form of the stochastic 
and dissipative interactions in a dissipative particle dynamics (DPD) model from projected 
phase space trajectories. The DPD model is viewed as a coai'se graining of a detailed 
dynamics that displays a clear time scale separation. Based on the Mori-Zwanzig projection 
operator method we deiive a consistency equation for the stochastic interaction in DPD. The 
consistency equation can be solved by an iterative boot strapping procedure. Combined 
with standard techniques for estimating the conservative interaction, our method makes it 
possible to reconstruct all the forces in a coarse grained DPD model. We demonstrate how the 
method works by recreating the interactions in a DPD model from its phase space trajectory. 
Furthermore, we discuss how our method can be used in realistic systems with finite time scale 
separation. 
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1. Introduction 

The molecular dynamics (MD) method is widely used to predict and analyze systems 
described by atomic or molecular models. Since only a small volume can be simulated, it 
is necessary to model how the simulated region interacts with the surroundings to bring the 
system into equilibrium. A device, or in the context of simulations a procedure, that brings 
the system to equilibrium is called a thermostat. Depending on the experimental setup the 
appropriate statistical model of the equilibrated system, the ensemble, is characterized either 
by constant energy and volume (NVE — MD without thermostat), by constant temperature 
and volume (NVT-e.g. the Nose-Hoover thermostat [1,2]), or by constant temperature and 
pressure (NPT-e.g. the Andersen [3] or Parrinello-Rahman [4] thermostat). Apart from 
attempting to mimic a specific experimental setup, the choice of thermostat might be guided 
for example by how easy it is to derive different thermodynamic properties. 

We can distinguish properties of the system that depend only on the equilibrium 
distribution of particle positions and velocities (e.g. the temperature, pressure, radial 
distribution function and the structure function) from transport properties (e.g. diffusion rate 
and viscosity). The former quantities depend only on the relative frequency of different states 
in the distribution, and because the dynamics of molecular systems is usually assumed to 
be ergodic, thermodynamic properties can be calculated from time averages over a single 
simulation run, rather than averaged over an explicit ensemble. 

Transport properties, contrary to equihbrium properties, are sensitive to the order 
in which states occur, i.e., to the temporal correlations. This follows from the Green- 
Kubo relations [5, 6], which give the transport properties in terms of autocorrelations of 
velocities and forces. The thermostat determines how the MD system approaches equilibrium, 
and its influence on the trajectories in turn alters the autocorrelations (compared to the 
constant energy dynamics) and therefore the transport properties. In addition to altering 
transport properties, most thermostats break fundamental symmetries of the systems, such as 
conservation of linear and angular momentum, since they are typically not Galilean invariant. 
This leads for instance to a damping of hydrodynamic modes [7]. In general, any effect of 
these thermostats on the statistics is therefore an artifact. In order to minimize these effects the 
strength of the thermostats are typically chosen small enough that the artifacts can be ignored, 
and the system size is chosen as large as possible. 

The effective dynamics of a coarse-grained molecular system has yet another source 
of dissipative interactions, compared to a system with only atomistic interactions. If the 
underlying system have a pronounced time scale separation, the thermostat can naturally 
appear as a result of projecting away the fast degrees of freedom [8,9]. MD itself may 
also be viewed as a coarse-grained representation of underlying quantum mechanics, based 
on the Born-Oppenheimer approximation. In standard MD the influence of the electron 
structure is usually neglected (except in ab initio MD), and the potentials are fitted to empirical 
data [10]. It is an open question whether a more formal coarse-graining would lead to a natural 
thermostat for MD systems. 

The theoretical framework underpinning this view of the origin of dissipative forces is 
Mori-Zwanzig theory of projection operators [8,9, 11-13]. Briefly, the theory explains how 
the dynamics of a microscopic system can be mapped to a coarse-grained or mesoscopic 
level by using projection operators. Depending on the projection, time-scale separation and 
the degree of exchange of energy between the fast and slow degrees of freedom, the effect 
of the fast degrees of freedom can either be eliminated due to averaging, resulting in a 
deterministic dynamics for the slow degrees of freedom, or it can result in Markovian white 
noise and dissipation, leading instead to a stochastic dynamics described by a Fokker-Planck 
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equation [14]. Used in equivalent Langevin-type equations, the drift and diffusion term in the 
Fokker-Planck equation naturally defines the thermostat for the coarse-grained system. 

An example of Langevin-type dynamics that can be coupled to the Mori-Zwanzig theory 
is the simulation technique dissipative particle dynamics (DPD). DPD is a particle based 
coarse grained model with pairwise central force interactions. The interactions have both a 
conservative part and a part given by noise and dissipation. By construction, the DPD model 
is Galilean invariant and can therefore be used to simulate nontrivial hydrodynamics. It was 
first suggested as a simulation tool for complex fluids [15], using soft conservative potentials. 
This has been the major conception of DPD, but it has also been shown that the method is 
suitable as an alternative thermostat for MD simulations [7]. 

An important conceptual shift is introduced if the dissipative and stochastic forces in 
DPD do not stem from interactions with the surroundings of the system but from interactions 
between the coarse-grained degrees of freedom and the degrees of freedom not explicitly 
modelled. In this situation, the thermostat should no longer be viewed as merely a means to 
make the system approach equilibrium, but as an integral part of the dynamics, where both 
the functional form and the strength of the thermostat are defined by the projection from the 
microscopic to the mesoscopic system. There exists formal coarse graining schemes resulting 
in DPD like mesoscopic dynamics [16, 17]. The main idea in this paper is to introduce a 
practical method for determining the functional form of the stochastic interactions in the DPD 
model from the phase space trajectories of the coarse grained system. We apply Mori-Zwanzig 
theory to derive a consistency expression that the stochastic interaction in the DPD model must 
fulfill. This result is used to derive a practical boot strapping method that can be used with 
simulation data to obtain a realistic estimate of the full functional forms of the effective coarse- 
grained interactions. In order to demonstrate the method and test its consistency, we apply it 
to phase space trajectories from a DPD simulation with known conservative, dissipative and 
stochastic forces. 

2. Theoretical analysis 

This section describes how the DPD model can be viewed as the effective dynamics resulting 
from a projection of an underlying atomistic dynamics. We begin with a general review of the 
projection operator method. We then discuss how DPD can be used as a specific ansatz for the 
effective dynamics resulting from center of mass projections of atomistic systems. We discuss 
the equilibrium and transport properties separately; First, for given stochastic forces, we show 
how the dissipative force must be chosen to maintain the equilibrium distribution, and how this 
relation leads to dissipative forces that respect the fundamental symmetries of the underlying 
dynamics. Second, we show that the combined effect of the stochastic and dissipative forces 
is to drive the system to thermal equilibrium, and that the radial dependence of the stochastic 
forces determines the rate of convergence to thermal equilibrium. In summary, this shows that 
the dissipative and stochastic parts of the effective coarse-grained dynamics are not arbitrary 
but are determined by the underlying dynamics through the choice of projection. 

2.1. Projection operators 
Consider a dynamical system 



X = f{x,y), 



(1) 
(2) 
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consisting of fast (y) and slow (x) degrees of freedom, and a time scale separation indicated 
by the parameter e ^ 1. The corresponding Liouville operator spUts into a fast and a slow 
part: 

- ^ = Y.i^J^^^^y)^\Y.^9^{x,y). (3) 



slow fast 
From the time evolution in density space, dtp = —Cp, an effective Fokker-Planck equation for 
the slow degrees of freedom can be derived. Following Just et al. [11], consider the adiabatic 
distribution pad{y\x) as the stationary distribution of y when x is considered fixed (changing 
adiabatically slowly). Under the assumption of ergodicity, an adiabatic average over the fast 
degrees of freedom conditioned on the slow degrees of freedom can be defined for an arbitrary 
function h, 

{h)ad{x)= [dy h{x,y)pad{y\x) = lim ^ [ dTh{x,T][T,y;x]), (4) 

J T~*oo 1 Jq 

where -qlr, y; x] is the trajectory of the differential equation r) = e^^g{x, ij) with 77(0) = y 
and X fixed. The ergodicity relation can be used in a practical situation to derive an adiabatic 
average from the detailed trajectory. The reduced phase space density can now be defined as 
p = J dypad{y\x), and the corresponding Fokker-Planck equation takes the form 



dpt _ ^ d n(l)^.^;T , ^ n(2) 
where the diffusion term is defined as 



1 i i' •' 



D^fix) = / dr {SfMx, 77[t, y; x])SFfjix, y)\^, (6) 
Jo 

and we use the notation 

hft{x, y) = Mx, y) - {S^)ad{x) (7) 

as an abbreviation for the fluctuations around the adiabatic equilibrium. The diffusion 
coefficient is given by the auto-correlation of the fluctuations in the fast degrees of freedom 
around their adiabatic stationary mean value. This is the relation that we will use to derive the 
functional form of the noise in DPD. At equilibrium, the drift term Z)'^^' is derived from the 
diffusion term using the fluctuation-dissipation theorem. 

In the form presented here, the Fokker-Planck equation represents the global dynamics 
on the phase space. There is no assumption that the system should have the structure of a 
mechanical system consisting of particles with pairwise interactions. For the methodology to 
be useful it is necessary to adopt it to a situation where an effective particle dynamics can be 
derived, e.g. a model with pairwise additive interactions like the DPD model. In Section|3]we 
continue the discussion on projection operators by showing how to apply the theory to derive 
a DPD dynamics from a more detailed, typically deterministic, simulation of a many-particle 
system. 
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2.2. The DPD model 

In a many-particle simulation a natural choice of coarse graining is to group particles in an 
underlying description of the system, e.g. an atomistic model, into single spherical particles 
(beads) positioned at the center of mass of the underlying particles. Furthermore, since forces 
between particles in the underlying system typically are pairwise and opposite, so that the 
system obeys Newton's third law, we assume that this also holds for the effective forces on 
the coarse-grained particles. This guarantees that the coarse-graining procedure does not 
break the conservation of linear and angular momentum, which in turn guarantees that the 
proper hydrodynamical behavior will be preserved in the approximation. The forces can be 
divided into three categories: conservative and dissipative deterministic forces, and stochastic 
forces. The conservative forces stem directly from the conservative interactions between the 
microscopic particles in one bead with the particles in another bead. The stochastic forces 
is the result of how fast chaotic degrees of freedom of the particles in each bead fluctuate 
around the motion of the center of mass and give rise to rapidly fluctuating forces. Finally, 
the dissipative forces represent the combined effect of the fast degrees of freedom on the slow 
degrees of freedom. With particles positioned at r^, with velocities and momenta p^, the 
equations of motion for a DPD model can be written as a system of Langevin equations 

p,=5^[Fg+Fg+F|.], (8) 

where F^, F^ and F^- are the conservative, dissipative, and stochastic forces between 
particles i and j. In DPD, the stochastic force between particles i and j take the form 

F,g = v/2fc^c^(ry)Gjey, (9) 

where j is the distance between particles i and j, is the unit vector pointing from j to i, 
fcs is Boltzmann's constant, and T is the temperature in Kelvin. The scalar function oj{rij) 
describes how the stochastic force depends on the distance between the particles, and Qj is 
interpreted as a symmetric Gaussian white noise term with mean zero and covariance 

{Cr3it)Qr{t')) = iS.u'S,j. + S,,.d,,,)d{t - t'), (10) 

where Sij and 5{t) are the Kronecker and Dirac delta functions, respectively. This structure 
of the covariance matrix makes sure that the stochastic forces between any pair of beads are 
central forces with equal magnitude, thus preserving the linear and angular momentum of the 
system. 

In order to illustrate how the DPD dynamics acts as a thermostat, and to further 
emphasize that the dissipative and stochastic parts of the dynamics are not arbitrary but 
are determined by the underlying dynamics through the choice of projection, we first derive 
how the damping force must be chosen in order to maintain the proper thermal equilibrium 
distribution of velocities and positions for a general Hamiltonian system, and then show 
how the radial dependence of the stochastic forces determine the rate of convergence to the 
equilibrium distribution, starting from an arbitrary distribution over phase space. 

2.3. Equilibrium dynamics 

At thermal equilibrium the system is distributed according to the canonical ensemble 

Pe,(x,p) = Z-lc-«(-P)/^-^, (11) 

where Z is the normalization term for the distribution. Since the equilibrium distribution 
depends only on the temperature and the Hamiltonian of the system, the conservative 
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forces uniquely determine the equilibrium distribution. The converse is also true; when the 
conservative forces depend only on the distance between particles, as in the DPD dynamics, 
it is possible to use an iterative approach to uniquely determine the conservative forces from 
the equilibrium radial distribution [18-23] (alternatively, direct time averaging over the fast 
degrees of freedom can be used, see e.g. Refs. [24, 25]). The methods for reconstructing 
effective potentials from the RDF rely on the result by Henderson [26], that two pairwise 
potentials resulting in the same RDF cannot differ by more than an additive constant. The 
importance of this theorem lies in the one-to-one correspondence between pairwise central 
force and the radial distribution function. 

The question is now: What is required of the forces to maintain the equilibrium 
distribution? The equilibrium ensemble is automatically invariant under the conservative parts 
of the dynamics, since i7 is a constant of motion for Hamiltonian dynamics (this is true for 
any ensemble where the probability of finding the system in a given micro-state depends 
only on the energy). The stochastic and dissipative forces generally change the equilibrium 
distribution, except if we choose the dissipative force such that the combined contributions 
from the dissipative and stochastic forces cancel when acting on the equilibrium distribution 
(the fluctuation-dissipation relation). Writing down the Fokker-Plank equation that describes 
the time evolution of the distribution over the state space in the DPD equations of motion, and 
requiring that the equilibrium distribution is a stationary point of the dynamics, leads to [27] 

= £^,--ff(x,p)/fcBT 



^ Vp, ■ [ - Ff + i ^ 2fcBrA,,(x)Vp^ 
V 



E 



Pi 



2 

J 
j 



-ff(x,p)/fcBT 



where C is the Fokker-Planck operator of Equation ([8]), and Ay is a 3 x 3 matrix given by the 
covariance of the total forces on particles i and j. The equilibrium Fokker-Planck equation 
( fT2] ) is commonly referred to as a fluctuation-dissipation relation. Since it must hold for all 
points (x, p) in phase-space, the only possible solution for the dissipative force is 

Ff = -^A,,(x)Vp^iJ. (13) 

3 

We briefly comment on the role of the coarse-graining projection in determining the 
dissipative forces. In DPD it is usually assumed that the projection is from a set of underlying 
particles to their center of mass. If this is not the case, but a more general projection is used, 
the equilibrium distribution of the projected dynamics is not necessarily the Gibbs distribution 
with a standard quadratic kinetic term. Especially, it may mean that the momentum-dependent 
part of the distribution depends also on space (i.e. that it is not possible to split the distribution 
into a momentum term cxp(— ^ • p^/2m,ifcsT) and a term that depends only on the positions 
of the DPD particles). This, in turn, means that the damping force may have a different 
shape than it has now (especially it may not be simply linear in the velocity). However, 
Equation ( fT3] l is still valid as long as there exist an energy function H such that the Gibbs 
distribution describes the equilibrium. 

For the DPD model the force covariance is given by 

(-[j^(rjj) ey (g) By when i ^ j 

uj^{nk)e,k<E)e,k wheni^j, '^^^^ 
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where ® denotes an outer product, and w(r) determines the radial dependence of the 
stochastic force. Inserting Equation ( fT4l ) into Equation ( fT3] l we can write the dissipative force 
on a particle as a sum of pair-wise dissipative forces: 

Ff - E = -Y.^\n3)^'^j ■ (v. - V,) e,,. (15) 

Equation ( fTsT i was first derived by [28]. Note that since Ff depends only on the velocity 
differences between interacting particles, it is manifestly Galilean invariant, and it is clear 
from the derivation above that this is a direct consequence of the covariance property of the 
stochastic forces [c.f. Equation (fTOlll. which in turn stems from the assumption that Newton's 
third law applies. 



2.4. Global convergence to equilibrium 

We have seen that the conservative part of the dynamics is determined by the equilibrium 
distribution, and the dissipative part of the dynamics is determined by the dependence of 
the Hamiltonian on the momentum of the beads in combination with the structure of the 
stochastic forces. This leaves only the stochastic forces to be determined in order to have 
a complete description of the DPD dynamics. The equilibrium distribution gives no hint 
here, since for any choice of stochastic force the dissipative force will maintain the Gibbs 
distribution. Instead, the choice of stochastic force will determine how the system approaches 
equilibrium. In order to better understand the effect of the stochastic forces on the path to 
thermal equilibrium, it is illuminating to study the time-evolution of the Kullback distance 
[29] from the non-equilibrium ensemble p(x,p) to the equilibrium distribution pe^lx, p), 
given by 

K{t) = / dxdpp(x,p) In^^^. (16) 

J Peq(x,p) 

The Kullback distance is non-negative for all ensemble distributions, and is zero if and 
only if the distribution is identical to the equilibrium distribution. With strictly Hamiltonian 
dynamics, K{t) is constant in time. Intuitively, this is because the internal energy of the 
system needs to change in order for the ensemble to approach the equilibrium distribution, 
but the Hamiltonian conserves the energy. Using the full Fokker-Planck equation of the DPD 
dynamics, we can calculate the rate of change of K{i): 

dtK{t) = -keT /dxdppE f ^P. log — ) f^P. log — 

■ /"dxdpp^ w(ry)e.y(Vp. ~ VpJ log — 

which is negative, except when the system is at equilibrium. Note that the change in K{t) 
does not depend directly on the conservative forces, but on the structure of the stochastic 
forces, uj{r). For any choice of uj{r), and from any initial distribution p over the phase space 
for the particle system, K{t) continues to decrease until p ~ p^q, where K = 0. 

Consider an initial distribution p(x, p) such that the system is in equilibrium with 
respect to the momentum space, but not with respect to position space (i.e. the momentum 
dependence of p and is the same, but the position dependence is different). In this case, 
we see that dtK = despite the system being out of equilibrium. This apparent paradox 
is resolved by calculating to the second derivative of K, to see that K is still a decreasing 
function of time {dfK < 0); such distributions correspond to saddle points in the dynamics. 



knT 



eg , 
2 

(17) 
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From this derivation it is apparent that the dissipative and stochastic forces act as a 
thermostat to bring the system to the equihbrium distribution, and that the rate at which this 
occurs, and by which path it occurs, is determined by the structure of the stochastic force (in 
combination with the conservative force). Thus, if we want the transport properties of our 
DPD system to match those of the underlying system, it is necessary to find the correct choice 
of stochastic force. In the next section, we describe a practical method for estimating the 
stochastic force function a;(r) from observed trajectories of the system. 

3. Estimating the stochastic force 

Coarse grained models of molecular systems are assumed to represent the projected dynamics 
of an underlying more detailed model. In the cases when the detailed dynamics can be 
simulated e.g. by the MD method, example trajectories of the projected system can be 
extracted by applying an explicit projection to the trajectories from the MD system. The 
detailed simulations of smaller systems over shorter time scales can then be used to calibrate, 
or as we do in this paper derive, the effective interactions in the coarse grained model. In 
this section we show that if the resulting reduced model fulfill the DPD ansatz (as formulated 
in Section \2.2\ . the observed trajectories contain enough information to estimate both the 
deterministic and stochastic forces in the DPD equations. As mentioned in Section |23] it 
is well established how the conservative force is determined from the equilibrium radial 
distribution function. In this section we propose a method for estimating w(r), determining 
both the dissipative and stochastic forces. 

3.1. Relation between DPD and the Mori-Zwanzig therory 

In order to relate the stochastic force in the DPD equations of motion to the observed 
behaviour of the particles, consider a pair of beads, i and j, a distance r from each other 
at time zero. As a first approach (and a naive one as we will explain), we expand the particle 
positions and momenta at a short time t to obtain 



where denotes a conditional ensemble average over all conformations where = r, and 
^tP = P(''') — p(0). For small r, only the leading term will be important, and can be used to 
estimate w^(r). 

Equation ( fTSl l relies explicitly on the infinite time-scale separation between the stochastic 
and deterministic forces. Such a separation can never hold in a natural system, since the coarse 
grained dynamics is smooth on the time scale of the underlying dynamics; this means that in 
the time region where Equation ( fTSb can be used to estimate cj(r), the DPD ansatz is not valid. 

A better approach follows from the Mori-Zwanzig coarse graining scheme laid out in 
Section im where the structure of the stochastic force was formulated in terms of a Green- 
Kubo relation (|6]l. To make contact between the general Mori-Zwanzig theory and the DPD 
model, the global dynamics is split into pairwise interactions, and the relative distances 
between pairs of particles and their relative velocities are assumed to be slow variables. Then, 
by comparing the stochastic term in the DPD model with ^ we find 



fcs-t Jo 

where SpYi is the difference between the projected total force on the beads and the adiabatic 
average, i.e.. Equation (|7]i adapted to the DPD ansatz. It is important to note that the adiabatic 



(18) 
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average is defined by all the deterministic forces on the reduced level, both conservative and 
dissipative. 

Rather than attempting to calculate the integral in Equation ( fT9l ) directly, we want to show 
how this can be estimated from observed trajectories of the coarse grained system. Using the 
notation 

^t(<5FP0 = / dt<5j.F,(t), (20) 
we find the identity 

2^ dt{5FV^{t)■5F¥,{Q))^ = ^{5t{5FP^)■5t{5FVJ)),■ (21) 

Thus, we can express Lo^{r) in terms of the asymptotic slope of the momentum change 
covariance: 

u'^ir) = lim -TTT^^ {5t{5FVi) ■ ^tiSpPj)), ■ (22) 

In practice it is not possible to take the limit to infinity, due to the so-called plateau 
problem [30]. This arises because for large values of r, the slope of {St{SpPi) ■ St{SpPj)) . 
must necessarily vanish, since eventually all beads separate and move independently. In 
addition, the measurements are conditioned on a distance r, rendering the right side of 
Equation (|22] | meaningless in the limit of large r since the constant r cannot be defined. 
If the fast and the slow parts of the dynamics are sufficiently well separated, however, 
there is a region where r is large enough that the fast degrees are in an approximate local 
(thermal) equilibrium determined by the slow degrees of freedom, but small enough that 
the slow degrees of freedom do not have time to change significantly. In this region, 
(StiSpPi) ■ StiSpPj)) is approximately linear, and the slope can be used to estimate the 
stochastic forces. If the conditioning on r does not hold exact, this can result in a small 
perturbation of {St{SFPi) ■ StiSFPj))^.- How to deal with this in practice will be explained in 
the next section. 



3.2. Boot strapping method 

The operator 5 f introduces a dependence of the dissipative force for the right side of Equation 
(l22b . The dissipative force in turn depends on Lo{r) (see Equation ( fTsT i) and Equation ( |22] | is 
therefore not closed. This can be resolved by a "boot strapping" approach. 

To estimate ijj{r) we assume access to a set of coarse grained trajectories, with the same 
time resolution as the underlying dynamics. uj{r) is found by solving Equation ( |22] | iteratively, 
with e.g. Lu{r) — Q in the first iteration. In the calculations, r is typically chosen to be much 
larger than the time steps in the underlying dynamics. 

The iteration procedure starts with the calculation of {5t{5FPi) ■ ^tiSpPj)) as a 
function of r for each value of r. To obtain an estimate of uj{r), the time region where 
the DPD ansatz is expected to be valid must be identified. This can be done by visual 
inspection, as illustrated in Figure [T] For small values of t (left section in Figure [U, the 
coarse grained dynamics follows the underlying dynamics smoothly and cannot be expressed 
in terms of DPD. Thereafter, the time region of interest follows (mid section in Figure[T]i. This 
region should be approximately linear in r, but this might not hold for two reasons. First, the 
measurements are conditioned on r being constant, but r is actually changing (slowly) with r. 
Second, unless the boot strapping procedure has converged, the dissipative force is not correct. 
Both these factors will affect the shape of the curves in Figure [T] To compensate for this, we 
fit the curves to a second order polynomial in the selected time region. Under the assumption 
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Figure 1. The figure shows a sketch of typical measurements of the momentum change 
covariance. Each curve (circles) con'esponds to a given value of r and for clarity only a few 
values of r are shown. In the left section, the dynamics of the coarse grained system is on the 
same time scale as the underlying dynamics and is therefore smooth. The mid section shows 
the region where the DPD ansatz is supposed to hold. This region should be approximately 
hnear, but can have higher order terms, as explained in the text. A second order polynomial fit 
is used in this region (solid lines). An estimate of u)'^ (r) is obtained from the coefficients of 
the hnear terms. 



that the second order terms are mainly a result of the conditioning on r (and the dissipative 
force when this has not yet converged), the coefficients of the linear terms give for each value 
of r the estimate of uj"^ (r). This procedure is repeated until convergence. As a word of caution, 
we recommend not to calculate the direct numerical derivative of the curves, as Equation ( l22b 
suggests. The numerical differentiation introduces noise and requires significantly longer 
simulations to obtain good statistics. It is much better to first do a (local) fit of the curve to a 
low-order polynomial and then evaluate the derivative of the fitted curve [31]. 



4. Numerical verification 



In order to verify the applicability of the proposed boot strapping method, we test the method 
on a DPD system with known ^(r-). The DPD simulations were set up using the standard 
implementation from Groot and Warren [32], with a density of 4.0 and a time step of 0.005. 
This time step is rather small, but we assume here that for a real coarse grained system, we 
will have access to the microscopic dynamics, which evolves on a time scale smaller than that 
usually used in DPD. The conservative force was chosen as 

pC ^ _j " " ' when r^j < 1, ^^^^ 




otherwise. 



with a — 25, and uj{r) was chosen as. 



Lo{r -) = I ~ """^^ ^^^"^ """^ ^ ^' (24) 
I otherwise, 
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Figure 2. Results from applying the boot strapping method in Section 13.21 to phase space 
trajectories from a DPD simulation. Curve A shows a;(r) after one iteration, starting with 
the initial guess u){r) = 0. Curves B and C show u){r) after the second and fourth iteration, 
respectively. The dashed line gives the exact value of u>(r) as used in the DPD simulation. 



with a — 3.0. These are standard parameter values chosen for the DPD fluid to match 
the compressibility of water at room temperature [32] and has been used extensively in 
mesoscopic simulations of lipid membranes [33,34]. 

As explained in the previous section, we wish to demonstrate that it is possible to 
obtain good estimates of uj{r) from Equation (l22l i. not only in the r ^ limit which is 
attainable in DPD but generally unavailable for a coarse grained system, but also for larger 
values of r. By using the boot strapping procedure outlined above for simulation data in the 
region r G [0.1,0.25], we show in Figure |2] the sequence of resulting uj{r) curves for the 
first boot strapping steps. Convergence towards the correct functional form of uj{r) is fast. 
Within the first four iterations, the method has approached the correct w-values used in the 
DPD simulation for most r-values. The conservative forces prevent the beads from coming 
arbitrarily close. This leads to a lack of data for small values of r, and explains the poor 
performance of the method in this region. For these values of r, we recommend to either 
set oj{r) to zero, or to use the value w(r*), where r* is the smallest value of r with reliable 
statistics in the simulations. 

To emphasize the importance of removing the full deterministic force (i.e. both 
conservative and dissipative forces) from the total force when calculating (5t(5FPi)' we show 
in Figure [3] the slope of the momentum change covariance, as defined in Equation (l22l i. for 
two different cases. In the first case (solid lines), Sp^i is defined as 



JfF, = F, - Ff - 
and in the second case (dashed lines), 

JfF, = F, - Ff , 



(25) 



(26) 



Fp is the conservative DPD force, and Ff* 
0, both methods converge to the correct 



where F; is the total force acting on particle i 
the dissipative DPD force. In the limit of t 
value of u;(r), as can be seen from Figure |3] but with Sp^i defined by Equation ( |25] |, there 
exists a plateau of small r-values for which the force covariance is approximately constant 
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0.05 0.1 0.15 0.2 0.25 

t 



Figure 3. The slope of the momentum change covariance as a function of r for different values 
of r. The curves are calculated from phase space trajectories of two DPD simulations, using 
different definitions of SpYi. The solid and dashed Unes correspond to using Equation (25) 
and Equation )26t . respectively. The different curves in the figure correspond to the r-values 
0.245 (top), 0.365, 0.485, 0.605 and 0.725 (bottom). The conservative force was set to in 
the simulations to demonstrate as clearly as possible the effect of not removing the dissipative 
force when calculating Sti^FPi) ■ 

Measuring (r) anywhere in this region gives approximately the correct value, and as shown 
previously, using r-values as far out as r £ [0.1, 0.25] still reproduces uj{r) using the boot 
strapping procedure. 

5. Discussion 

We have discussed how to derive a dissipative particle dynamics from a detailed microscopic 
system, for example a molecular dynamics simulation. As a coarse grained model of 
a mechanical system, DPD has several advantages compared to many other models. 
Most importantly, by construction the DPD dynamics respects fundamental symmetries of 
the underlying dynamics; it is Galilean invariant, and therefore both linear and angular 
momentum are locally conserved by the interactions. 

In this paper, we have established two important properties of our method: First, that 
the framework is consistent; if the dynamics is on the DPD form, we can use the method to 
accurately reconstruct all terms in the equations of motion. Through our adaption of the Mori- 
Zwanzig projection operator methods, we argue that this provides a clear and quantifiable 
connection to the underlying degrees of freedom. 

Second, the method successfully reconstructs the dynamics without using the shortest 
time-scale of the particle trajectories. If the detailed model shows a strong time-scale 
separation it is possible to use Equation ( fTsT i directly to estimate the effective stochastic 
interactions on the coarse grained level. In many cases however - e.g. when DPD is used 
as a coarse-grained representation of a molecular system - the time scale separation is not 
very significant. It is therefore essential that our method is able to reconstruct the dynamics 
without using the correlations of the system at the shortest time scale, which we demonstrate 
in Figure [T] 
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The price we pay for not using the short-time properties is that we cannot use a direct 
method to measure the shape of a;(r), but are forced to use an iterative scheme. This is 
because u){r) is estimated from the stochastic force, and in order to extract this force from 
the dynamics we need to subtract influence of the deterministic (conservative and dissipative) 
forces from the particle trajectories, which then depend on the u){r) that we try to estimate in 
the first place. However, we have found that starting from the initial u){r) — 0, and using the 
estimation as a fixed-point scheme, leads to rapid convergence (Figure|2]i; the likely cause for 
this is that the dissipative component of the deterministic force is typically dominated by the 
conservative component, so that akeady the first iteration leads to an Lo{r) not far from the 
correct one. 

The main advantages of our method is perhaps that it gives the appropriate magnitude 
of dissipative and stochastic forces for the coarse-grained system to be consistent with the 
underlying dynamics; hence, if the coarse grained dynamics is averaging, so that fluctuations 
are not important, the resulting ijj{r) k, 0, and if the rapidly fluctuating degrees of freedom act 
as white noise on the coarse-grained dynamics, Lo{r) captures this effect. This is in contrast 
to most thermostats used in MD, where in principle the thermostat is used to stabiUze the 
dynamics and is not considered to be an integral part of the dynamics. 
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